Material-based subdomain hybrid cellular automata algorithm for material optimization of thin-walled frame structures

ABSTRACT

A material-based subdomain hybrid cellular automata algorithm for solving material optimization of thin-walled frame structures includes an outer loop and an inner loop. The outer loop is to define and update the target cost for the inner loop, and the inner loop is to adjust material using a PID control strategy according to the nominal flow stress of a current cell and the nominal flow stress of candidate materials, so that a current cost of the inner loop converges to the target cost. The material-based subdomain hybrid cellular automata algorithm can effectively solve a nonlinear dynamic response optimization problem in the discrete design space.

CROSS REFERENCE TO THE RELATED APPLICATIONS

This application is the national phase entry of International Application No. PCT/CN2021/082325, filed on Mar. 23, 2021, which is based upon and claims priority to Chinese Patent Application No. 202110238581.2, filed on Mar. 4, 2021, the entire contents of which are incorporated herein by reference.

TECHNICAL FIELD

The present invention relates to the technical field of crashworthiness optimization of thin-walled frame structures, and specifically to a material-based subdomain hybrid cellular automata algorithm for solving material optimization of thin-walled frame structures.

BACKGROUND

Bodies of car, rail transit and engineering machinery are typical space frame structures assembled by a variety of thin-walled structures of different materials. To use multi-material thin-walled frame structures is an effective way and an inevitable trend to optimize the safety and cost requirements of automotive, rail transit and engineering machinery. To select the best material and to do optimization design of the multi-material thin-walled frame structures can not only improve the mechanical properties (such as crashworthiness) of thin-walled frame structures, but also reduce their costs. The multi-material matching optimization of thin-walled frame structures with considering crashworthiness and cost is a typical nonlinear dynamic response optimization problem including many discrete variables. Due to the high level of non-linearity in the output response of crash simulation and the presence of numerical noise and oscillation, the gradient-based optimization algorithm is difficult to effectively solve the crashworthiness optimization problem of thin-walled frame structures. Evolutionary algorithms usually need to perform thousands of finite element analyses (FEAs), which leads to a very long optimization time. The optimization algorithms based on the surrogate model is the main way to solve the problems mentioned above. However, when the number of design variables is large (such as more than 30 or even more), the optimization efficiency of the optimization algorithm based on the surrogate model will also be greatly reduced.

Hybrid Cellular Automata (HCA) method is a non-gradient heuristic algorithm that can solve nonlinear dynamic response optimization problems including many discrete variables (such as material density or thickness). However, the existing methods are mainly based on the idea of uniform distribution of internal energy density (IED) to update material density or thickness. It is difficult to solve the multi-material optimization problem of thin-walled frame structures with multiple performance constraint functions, and easy to fall into local optimal solution. Presently, few efficient algorithms can be employed to solve the nonlinear dynamic response optimization problem with many material variables in the discrete design space.

SUMMARY

In view of the deficiencies in the prior art, the present invention provides a material-based subdomain hybrid cellular automata algorithm for solving material optimization of thin-walled frame structures which can efficiently solve a nonlinear dynamic response optimization problem with many material variables.

The present invention achieves the technical object mentioned above by the following technical means.

A material-based subdomain hybrid cellular automata(M-SHCA) algorithm for solving material optimization of thin-walled frame structures includes the following steps:

-   -   S1. establishing an initial crash finite element model,         constructing a subdomain cellular automata model, defining         material variables and field variables of thin-walled frame         structures, and employing the initial crash finite element model         for material and cost optimization;     -   S2. executing an outer loop: calculating a cell IED and a         constraint value at a current design point by FEA, and updating         the target cost by the penalty function method according to the         extent of the current design point violating the constraint         boundary;     -   S3. executing an inner loop with the following steps:     -   S3.1. constructing a step IED target (SIED*) function and         updating the target IED;     -   S3.2. updating cell material using material updating rule based         on a PID control strategy;

Specifically: defining a candidate material library and a nominal flow stress of each material, updating the nominal flow stress of a current cell, comparing a nominal flow stress with the true flow stress of each material in the candidate material library, selecting a candidate material closest to the nominal flow stress as a selected material of current cell, and replacing material parameters of current cell with the mechanical parameters of the selected material;

-   -   S3.3. processing to S4 if an inner loop is convergent, otherwise         returning to S3.1;     -   S4. writing out an optimization result if the global convergence         conditions in the outer loop are satisfied, otherwise returning         to S2 for updating cell material in the inner loop.         Further, a candidate material library is defined as follows:

$\begin{matrix} {{{Mat} = \left\{ {{{Mat}(1)},L,{{Mat}(s)},L,{{Mat}(1)}} \right\}},{1 \leq s \leq l}} \\ {= \begin{matrix} \left\{ {\left( {\rho_{1},E_{1},\sigma_{y1},\sigma_{u1},\sigma_{f1},L} \right),L,\left( {\rho_{s},E_{s},\sigma_{ys},} \right.} \right. \\ \left. {{\left. {\sigma_{us},\sigma_{fs},L} \right)L},\left( {\rho_{l},E_{l},\sigma_{yl},\sigma_{ul},\sigma_{fl},L} \right)} \right\} \end{matrix}} \end{matrix}$

where, ρ_(s) is a density of the sth material in the candidate material library; E_(s) is an elastic modulus of the sth material in the candidate material library; σ_(ys) is a yield strength of the sth material in the candidate material library;

$\sigma_{fs} = \sqrt{\frac{\sigma_{ys}\sigma_{us}}{1 + n}}$

is a flow stress of the sth material in the candidate material library; σ_(us) is a tensile strength of the sth material in the candidate material library; :1≥2 is the number of materials in the candidate library. Furthermore, a nominal flow stress is a non-physical parameter, which is a positive real number. Further, a nominal flow stress of each cell is updated as follows:

$\sigma_{\Omega_{i,j}}^{\prime({{h + 1},k})} = \left\{ \begin{matrix} {\sigma_{\Omega_{i,j}}^{\min},} & {{\sigma_{\Omega_{i,j}}^{\prime({h,k})} + {\Delta\sigma}_{\Omega_{i,j}}^{\prime({h,k})}} < \sigma_{\Omega_{i,j}}^{\min}} \\ {{\sigma_{\Omega_{i,j}}^{\prime({h,k})} + {\Delta\sigma}_{\Omega_{i,j}}^{\prime({h,k})}},} & {\sigma_{\Omega_{i,j}}^{\min} \leq {\sigma_{\Omega_{i,j}}^{\prime({h,k})} + {\Delta\sigma}_{\Omega_{i,j}}^{\prime({h,k})}} \leq \sigma_{\Omega_{i,j}}^{\max}} \\ {\sigma_{\Omega_{i,j}}^{\max},} & {{\sigma_{\Omega_{i,j}}^{\prime({h,k})} + {\Delta\sigma}_{\Omega_{i,j}}^{\prime({h,k})}} > \sigma_{\Omega_{i,j}}^{\max}} \end{matrix} \right.$

where, σ_(Ω) _(i,j) ′^((h,k)) is a nominal flow stress of the jth cell in the subdomain Ω₁ in the hth inner loop and the kth outer loop; σ_(Ω) _(i,j) ′^((h+1,k)) is a nominal flow stress of the jth cell in the subdomain Ω₁ in the (h+1)th inner loop and the kth outer loop; σ_(Ω) _(i,j) ^(min) and σ_(Ω) _(i,j) ^(max) are a minimize value and a maximum value of actual flow stress of the jth cell in the subdomain Ω₁, respectively; Δσ_(Ω) _(i,j) ′^((h,k)) is a nominal flow stress variation of the jth cell in the subdomain Ω₁ in the hth inner loop and the kth outer loop:

Δσ_(Ω) _(i,j) ′^((h,k))=(σ_(Ω) _(i,j) ^(max)−σ_(Ω) _(i,j) ^(min))f(e _(Ω) _(i,j) ^((h,k)))

where, e_(Ω) _(i,j) ^((h,k)) denotes a difference between the current IED S_(Ω) _(i,j) ^((k)) and a target IED S_(m′)*^((h,k)) and a PID control function for updating the nominal flow stress f(e_(Ω) _(i,j) ^((h,k))) is given as follows:

${f\left( e_{\Omega_{i,j}}^{({h,k})} \right)} = {{K_{p}e_{\Omega_{i,j}}^{({h,k})}} + {K_{i}\left\lbrack {e_{\Omega_{i,j}}^{({h,k})} + {\sum\limits_{\tau = 1}^{k - 1}e_{\Omega_{i,j}}^{(\tau)}}} \right\rbrack} + {K_{d}\left\lbrack {e_{\Omega_{i,j}}^{({h,k})} - e_{\Omega_{i,j}}^{({k - 1})}} \right\rbrack}}$

where, K_(p) is a proportional control coefficient; K_(i) is an integral control coefficient; K_(d) is a differential control coefficient; e_(Ω) _(i,j) ^((r)) is a relative deviation item of the rth outer loop; e_(Ω) _(i,j) ^((k−1)) is a relative deviation item of the (k−1)th outer loop. Furthermore, the following equation is used to select a candidate material closest to the nominal flow stress as the selected material of the current cell by comparing the nominal flow stress with the true flow stress of each material in the candidate material library:

$\left\{ \begin{matrix} {{\sigma_{\Omega_{i,j}}^{\prime({{h + 1},k})} - \sigma_{fp}} = {\min\left( {{\sigma_{\Omega_{i,j}}^{\prime({{h + 1},k})} - \sigma_{f1}},L,{\sigma_{\Omega_{i,j}}^{\prime({{h + 1},k})} - \sigma_{fs}},L,{\sigma_{\Omega_{i,j}}^{\prime({{h + 1},k})} - \sigma_{fl}}} \right)}} \\ {\sigma_{\Omega_{i,j}}^{({{h + 1},k})} = \sigma_{fp}} \\ {{Mat}_{\Omega_{i,j}}^{({{h + 1},k})} = {{Mat}(p)}} \end{matrix} \right.$

where, p denotes a position of the selected material in the candidate material library, σ_(fp) is an actual flow stress of the selected material; σ_(Ω) _(i,j) ^((h+1,k)) is an actual flow stress of the jth cell in the subdomain Ω₁ in the (h+1)th inner loop and the kth outer loop; Mat(p) denotes a selected material; Mat_(Ω) _(i,j) ^((h+1,k)) denotes a selected material of the jth cell in the subdomain Ω₁ in the (h+1)th inner loop and the kth outer loop; σ_(s) is a flow stress of the sth material in the candidate material library. Further, the following equation is employed to replace material parameters of a current cell with mechanical parameters of selected material:

$\begin{pmatrix} \rho_{\Omega_{i,j}}^{({{h + 1},k})} \\ E_{\Omega_{i,j}}^{({{h + 1},k})} \\ \sigma_{y,\Omega_{i,j}}^{({{h + 1},k})} \\ \sigma_{u,\Omega_{i,j}}^{({{h + 1},k})} \\ \sigma_{f,\Omega_{i,j}}^{({{h + 1},k})} \\ L \end{pmatrix} = \begin{pmatrix} \rho_{p} \\ E_{p} \\ \sigma_{yp} \\ \sigma_{up} \\ \sigma_{fp} \\ L \end{pmatrix}$

where, ρ_(Ω) _(i,j) ^((h+1,k)) is a material density of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop; E_(Ω) _(i,j) ^((h+1,k)) is an elastic modulus of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop; σ_(y,Ω) _(i,j) ^((h+1,k)) is a yield tress of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop; σ_(u,Ω) _(i,j) ^((h+1,k)) is a tensile strength of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop; σ_(f,Ω) _(i,j) ^((h+1,k)) is a flow stress of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop; ρ_(p) is a material density of a selected material; E_(p) is an elastic modulus of the selected material; σ_(yp) is a yield tress of the selected material; σ_(up) is a tensile strength of the selected material; σ_(fp) is a flow stress of the selected material; σ_(up) is a tensile strength of selected material; σ_(fp) is a flow stress of selected material.

Furthermore, global convergence conditions comprise:

|C ^((h,k)) −C* ^((k))|<ε₁ or k ₁ ≥k _(1max)

where

$C^{({h,k})} = {\overset{l}{\sum\limits_{i = 1}}{\overset{N}{\underset{j=1}{\sum}}C_{\Omega_{i,j}}^{({h,k})}}}$

is a total cost in the hth inner loop and the kth outer loop, here C_(Ω) _(i,j) ^((h,k)) is a cell cost in the hth inner loop and the kth outer loop; C*^((k)) is a target cost defined in the kth outer loop; ε₁ is an inner loop convergence factor; k₁ is the number of iterations in the inner loop; k_(1max) is a maximum number of iterations in the inner loop; a cell cost of the selected material of the jth cell in the subdornain Ω₁ in the (h+1)th inner loop and the kth outer loop C_(Ω) _(i,j) ^((h+1,k)) is calculated as follows:

C _(Ω) _(i,j) ^((h+1,k))=ξ_(Ω) _(i,j) ^((h+1,k))ρ_(Ω) _(i,j) ^((h+1,k)) t _(Ω) _(i,j) A _(Ω) _(i,j)

where, ξ_(Ω) _(i,j) ^((h+1,k)) is a price of the selected material of the jth cell in the subdomain Ω₁ in the (h+1)th inner loop and the kth outer loop, ρ₁₀₆ _(i,j) ^((h+1,k)) is a density of the selected material of the jth cell in the subdomain Ω₁ in the (h+1)th inner loop and the kth outer loop, t_(Ω) _(i,j) is a thickness of the jth cell in the subdomain Ω₁, and A_(Ω) _(i,j) is an area of the jth cell in the subdomain Ω₁.

The present invention has the following beneficial effects:

-   -   (1) The inner loop introduces the nominal flow stress in the         present invention which converts discrete material variables         into continuous variables, and then the cell material updating         rule based on a PID control strategy is used to realize the         iterative updating of cell material and improve the robustness         of the M-SHCA algorithm.     -   (2) The M-SHCA algorithm can solve nonlinear dynamic response         optimization problems containing many material variables which         not only has high searching efficiency (less computational         resources and less finite element simulations), but also has         good global convergence accuracy.     -   (3) The present invention introduces the concept of subdomain         cellular automata model based on the topological structural         characteristics of thin-walled frame structure to solve the         nonlinear dynamic response optimization problems in a discrete         design space.     -   (4) A step target IED update rule is used in the present         invention to update the cell material in the inner loop, thereby         improving the global optimal solution searching ability of the         M-SHCA. algorithm.     -   (5) Gradient information is not necessary to be calculated         during the optimization process of the present invention, which         has a great advantage in solving complex nonlinear problems         where sensitivity information is hard to obtain.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flowchart of the M-SHCA algorithm for solving material optimization of thin-walled frame structure in the present invention;

FIG. 2 is a schematic diagram of the definition of a subdomain cellular automata (SCA) model in the present invention;

FIG. 3 is a schematic diagram of a horizontal internal energy density target (HIED*) function in the present invention;

FIGS. 4A-4B are schematic diagrams of a step target IED function in the present invention, in which, FIG. 4A is a schematic diagram of a step target IED function when H_(threshold)=0 and V_(threshold)=1 and FIG. 4B is a schematic diagram of a step target IED function when H_(threshold)=3 and V_(threshold)=1.1;

FIG. 5 is a schematic diagram of a full-vehicle side crash finite element model in the present invention;

FIG. 6 is a schematic diagram of a SCA model of a car body beam frame in the present invention;

FIG. 7 is effective stress-strain curves of candidate materials in the present invention;

FIG. 8 is measuring points at B-pillar, in which points A and B are respectively locations related to the loading positions to passenger's chest and pelvis;

FIGS. 9A-9C are curve diagrams of iteration processes of M-HCA#1 algorithm in the present invention, where FIG. 9A is a diagram illustrating total cost convergence curve, FIG. 9B is a diagram of an iteration process with a maximum intrusion amount, and FIG. 9C is a diagram of an iteration process with a maximum intrusion velocity;

FIGS. 10A-10C are curve diagrams of iteration processes of M-HCA#2 algorithm in the present invention, where FIG. 10A is a diagram illustrating total cost convergence curve, FIG. 10B is a diagram of an iteration process with a maximum intrusion amount and FIG. 10C is a diagram of an iteration process with a maximum intrusion velocity;

FIGS. 11A-11B are curve diagrams of iteration processes of a parallel EGO-PCEI algorithm in the present invention, where FIG. 11A is a diagram illustrating total cost convergence curve, FIG. 11B is a diagram of an iteration process with a. maximum intrusion amount, and FIG. 11C is a diagram of an iteration process with a maximum intrusion velocity;

FIG. 12 is a schematic diagram of comparing the computational resources and computational time of M-HCA#1, M-HCA#2 and parallel EGO-PCEI algorithms;

FIGS. 13A-13D are comparison diagrams of an intrusion amounts and an intrusion velocity of measuring points at B-pillar before and after optimization in the present invention, where FIG. 13A is a comparison diagram of B-pillar intrusion amount curves at chest location, FIG. 13B is a comparison diagram of B-pillar intrusion amount at pelvis location, FIG. 13C is a comparison diagram of B-pillar intrusion velocity at chest location, and FIG. 13D is a comparison diagram of B-pillar intrusion velocity at pelvis location;

FIGS. 14A-14B are comparison diagrams of deformation patterns of a car body before and after optimization in the present invention, where FIG. 14A is a diagram of a car body deformation pattern before optimization, and FIG. 14B is a diagram of a car body deformation pattern after optimization.

DETAILED DESCRIPTION OF THE EMBODIMENTS

The present invention is further illustrated below with reference to the accompanying drawings and specific embodiments, but the protection scope of the present invention is not limited thereto.

As shown in FIG. 1 , a material-based subdomain hybrid cellular automata(M-SHCA) algorithm for solving material optimization of thin-walled frame structures specifically includes the following steps:

-   -   S1: An initial crash finite element model of a thin-walled frame         structures is established for the material and cost optimization         of thin-walled frame structure;

A finite element preprocessor software (such as HyperMesh or LS-Prepost) is used to discrete a full-vehicle geometric model into its finite element meshes, and then assign attributes, materials, boundary conditions and initial conditions for each part of the finite element meshes to complete the full-vehicle crash finite element model for the material and cost optimization of thin-walled frame structures.

-   -   S2: A subdomain cellular automata model of thin-walled frame         structures is established and its material variables and field         variables are defined;

The concept of “subdomain Cellular Automata (CA) model” is introduced in a discrete design space based on a conventional CA model, as shown in FIG. 2 . Assume that a global design space Ω is fragmented into/discrete subdomains Ω_(i) (i+1,2,L,l) and each subdomain Ω₁ is described by a cellular automata model CA₁₀₆ _(i) and a cell state α₁₀₆ _(i) by the following mathematical expression:

Ω(CA _(Ω),α_(Ω))=Ω₁(CA _(Ω) ₁ ,α_(Ω) ₁ )+Ω₂(CA _(Ω) ₂ ,α_(Ω) ₂ )+L+Ω _(i)(CA _(Ω) _(i) ,α_(Ω) _(i) )+L+Ω ₁(CA _(Ω) ₁ ,α_(Ω) ₁ )   (1)

where Ω₁ denotes an ith subdomain of the global design space Ω, α_(Ω) _(i) denotes a cell state of the ith subdomain, CA_(Ω) _(i) is the cellular automata model of the ith subdomain, which is consisted of a current cell and its neighboring cells thereof, and the cell type maybe one-dimensional cellular automata, two-dimensional cellular automata, and three-dimensional cellular automata.

As shown in FIG. 2 , a jth cell state a α₁₀₆ _(i,j) in the ith subdomain Ω₁ is described by the following equation:

$\begin{matrix} {\alpha_{\Omega_{i,j}} = \begin{Bmatrix} {Mat}_{\Omega_{i,j}} \\ S_{\Omega_{i,j}} \end{Bmatrix}} & (2) \end{matrix}$

in which, Mat_(Ω) _(i,j) is a design variable (such as, part material) of a jth cell in the ith subdomain Ω_(i); S_(Ω) _(i,j) is a field variable (such as, part IED) of a jth cell in the ith subdomain Ω₁, calculated as follows:

$\begin{matrix} {S_{\Omega_{i,j}} = {\frac{1}{{\hat{N}}_{\Omega_{i,j}} + 1}{\sum}_{n \in N_{\Omega_{i,j}}}\frac{U_{\Omega_{i,n}}}{t_{\Omega_{i,n}} \cdot A_{\Omega_{i,n}}}}} & (3) \end{matrix}$

Wherein, U_(Ω) _(i,n) is an internal energy of the nth cell in the ith subdomain Ω₁, t_(Ω) _(i,n) and A_(Ω) _(i,n) are respectively a thickness and a surface area of the nth cell in the ith subdomain Ω₁, N_(Ω) _(i,j) is a set of the neighboring cells of the jth cell in the ith subdomain, Ω₁, S₁₀₆ _(i,j) is a local IED of the jth cell in the ith subdomain Ω₁, namely, the field variable, and {circumflex over (N)}_(Ω) _(i,j) denotes the number of cells in the collection.

-   -   S3: The process enters an outer loop, the finite element         analysis software, such as LS-DYNA, Radioss, Abaqus or Ansys, is         invoked to conduct a nonlinear dynamic crash simulation         analysis, a cell IED and a constraint function value at a         current design point are accessed from the simulation result,         and a target cost is updated by a penalty function method         according to an extent to which the current design point         violates a constraint boundary;

A main purpose of the outer loop is to carry out finite element simulation analysis, calculate an output response, and update an TED and a target cost:

-   -   S3.1: The IED updating rule

To improve the stability for updating the cell material and avoid the oscillation in the outer loop, an IED S_(Ω) _(i,j) ^((k)) is updated with the weighted IEDs of the last three iterations in the outer loop:

$\begin{matrix} \left\{ \begin{matrix} {{{\hat{S}}_{\Omega_{i,j}}^{(1)} = S_{\Omega_{i,j}}^{(1)}},} & {k = 1} \\ {{{\hat{S}}_{\Omega_{i,j}}^{(2)} = {{\frac{1}{2}S_{\Omega_{i,j}}^{(1)}} + {\frac{1}{2}S_{\Omega_{i,j}}^{(2)}}}},} & {k = 2} \\ {{{\hat{S}}_{\Omega_{i,j}}^{(k)} = {{\frac{1}{2}S_{\Omega_{i,j}}^{(k)}} + {\frac{1}{3}S_{\Omega_{i,j}}^{({k - 1})}} + {\frac{1}{6}S_{\Omega_{i,j}}^{({k - 2})}}}},} & {k \geq 3} \end{matrix} \right. & (4) \end{matrix}$

-   -   S3.2: A target cost updating rule

When multiple performance constraint functions exist in the outer loop, a target cost updating rule is proposed based on the penalty function method, in which the penalty value of a target cost ΔC*^((k)) is used to indicate an extent to a current design point violating the constraint boundary in the kth outer loop and then the target cost is updated in the kth outer loop. The penalty value of the target cost ΔC*^((k)) is calculated as follows:

$\begin{matrix} {{\Delta C^{*{(k)}}} = {\min\left( {{\max\left( {{K_{q}C^{*{(0)}}\delta^{(k)}},{{- \Delta}C}} \right)},{\Delta C}} \right)}} & (5) \end{matrix}$ $\begin{matrix} {\delta^{(k)} = {\sum\limits_{i = 1}^{n_{g}}\frac{O_{i}^{(k)} - O_{i}^{*}}{O_{i}^{*}}}} & (6) \end{matrix}$

in which, n_(g) is the number of constraint functions, K_(q) is a scale factor of the penalty value of the target cost, O_(i) ^((k)) is a response of the ith constraint function in the kth outer loop, O_(i)* is a specified boundary condition of the ith constraint function, δ^((k)) is a relative deviation between n_(g) constraint functions and a specified boundary condition, C*⁽⁰⁾ denotes an initial total cost of thin-walled frame structures, ΔC denotes a maximum penalty of the target cost.

Then a target cost C*^((k)) in the kth outer loop is updated by the following equation:

C* ^((k))=min(C* ^((k−l))+ΔC* ^((k)) ,C* ^((k′)))   (7)

in which, k′ denotes the position of the last feasible solution in the outer loop iterations. If there is no feasible solution in the outer loop iterations until current design point, k′ will vanish (k′=0).

To improve the convergence efficiency of the M-SHCA. algorithm, p_(f) is defined to indicate the iteration number for the consecutive infeasible solutions, of which an initial value is set to be zero; p_(f)* is defined to indicate the maximum iteration number for the consecutive infeasible solutions. If the current design point is a feasible solution during iterations, p_(f)=0; if the current design point is an infeasible solution, p_(f)=p_(f)+1. If the iteration number for the consecutive infeasible solutions is greater than the maximum iteration number for the consecutive infeasible solutions (p_(f)>p_(f)*), the M-SHCA algorithm will be convergent and the iterations of the M-SHCA algorithm will be terminated.

-   -   S4: The process enters an inner loop, and the following steps         are performed:     -   S4.1: A step target internal energy density function (step IED         target, SIED*) is constructed, and a target internal energy         density is updated;

As shown in FIG. 3 , a so-called horizontal IED target (HIED*) function is employed to update the cell material which will keep increasing to an upper limit for the cells in the high IED subdomains or decreasing to a lower limit for the cells in the low IED regions. However, the IEDs in the whole design space are still very difficult to converge to HIED*. So, we construct a step IED target (SIED*) function with the IED targets of different subdomains calculated adaptively by the cells' IED distribution. Specifically, according to the cells' IED distribution, multiple “step points” and “step intervals” are adaptively defined and a reasonable IED target is calculated for each “step interval”, in which the cell material is updated by the SIED* function, and the searching robustness of the global optimal solution is improved.

A side collision simulation of car body frames is employed to exemplify the construction and the updating rule of the SIED* function:

-   -   S4.1.1: Cell index definition

An index id of the cell Ω_(i,j) with its subscripts i and j is defined by equation (8), that is, id is a function about the subscripts i and j of the cell Ω_(i,j), and then S_(id(i,j)) ^((k))=S₁₀₆ _(i,j) ^((k)),

id(i,j)={circumflex over (N)}_(Ω) _(i−1) *(i−1)+j, (j∈[1,{circumflex over (N)}_(Ω) _(i) ], {circumflex over (N)}_(Ω) ₀ =0)   (8)

in which, {circumflex over (N)}_(Ω) _(i−1) is the number of cells in the (i−1)th subdomain. A specific case about the cell index definition is given in TABLE 1.

TABLE 1 Correspondence relationship between the cell Ω_(i, j) and its index Cell Ω_(1, 1) Ω_(1, 2) Ω_(1, 3) Ω_(1, 4) Ω_(2, 1) Ω_(2, 2) Ω_(3, 1) . . . Ω_(11, 1) Ω_(12, 1) Ω_(13, 1) Ω_(14, 1) Index 1 2 3 4 5 6 7 . . . 21 32 33 34

-   -   S4.1.2: All cells are traversed, and the differences between the         IED S_(id) ^((k)) of each cell and an average IED S ^((k)) of         all cells in the kth outer loop are calculated:

ΔS _(id) ^((k)) =S _(id) ^((k))−S ^((k))   (9)

where,

${\overset{\_}{S}}^{(k)} = {\frac{1}{\sum\limits_{i = 1}^{l}{\hat{N}}_{\Omega_{i}}}{\sum\limits_{i = 1}^{l}{\sum\limits_{j = 1}^{{\hat{N}}_{\Omega_{i}}}S_{\Omega_{i,j}}^{(k)}}}}$

is an average IED of all cells in the kth outer loop,

-   -   S4.1.3: Determination of “step points” and “step ranges”

All cells are traversed to judge whether equation (10) is satisfied. A subscript id of ΔS_(id) ^((k)) is defined as a “step point” and denoted as id_(i) if equation (10) is satisfied. The m “step points” determined by equation (10) can construct m+1 “step ranges” denoted as [id_(i−1), id_(i)], where

${i = 1},\ldots,{m + 1},{{id_{0}} = 1},{{id}_{m + 1}{\sum\limits_{i = 1}^{l}{{\hat{N}}_{\Omega_{i}}.}}}$ ΔS_(id) ^((k))*ΔS_(id+1) ^((k))<0   (10)

-   -   S4.1.4: Update of “step points” and “step ranges”

A width threshold of the “step range” is defined as H_(threshold). All “step ranges” are traversed. to judge whether equation (11) is satisfied. If equation (11) is satisfied (that is, the width of “step range” [id_(i−1), id_(i)] is very small), the “step points” are deleted and the “step ranges” are updated in the following manner: when i=1, a “step point” id₁ is deleted, the “step range” is updated from [id₀,id₁] to [id₀,id₂]; when i>1, a “step point” id_(i−1) is deleted, and the “step range” is updated from [id_(i−1),id_(i)] to [id_(i−2),id_(i)]. The original “step points” and “step ranges” are retained if equation (11) is not satisfied. If the number of the updated “step points” is m′, the number of the updated “step ranges” is m′+1.

id_(i+1)−id_(i)+1<H_(threshold)   (11)

-   -   S4.1.5: A step target IED function is constructed as follows:

$\begin{matrix} {S^{*{({h,k})}} = \left\{ \begin{matrix} {S_{1}^{*{({h,k})}},{1 \leq {id} \leq {id}_{1}}} \\ {S_{2}^{*{({h,k})}},{{id}_{1} \leq {id} \leq {id}_{2}}} \\ {\ldots} \\ {S_{i}^{*{({h,k})}},{{id}_{i - 1} \leq {id} \leq {id}_{i}}} \\ {\ldots} \\ {S_{m^{\prime}}^{*{({h,k})}},{{id}_{m^{\prime} - 1} \leq {id} \leq {id}_{m^{\prime}}}} \\ {S_{m^{\prime} + 1}^{*{({h,k})}},{{id}_{m^{\prime}} \leq {id} \leq {\sum\limits_{i = 1}^{l}{\hat{N}}_{\Omega_{i}}}}} \end{matrix} \right.} & (12) \end{matrix}$

where, S_(i)*^((h,k)) is a target IED in the “step range”[id_(i−1),id_(i)] in the kth outer loop and the hth inner loop.

-   -   S4.1.6: Update of the step target IED function

To achieve the specified target mass in the outer loop, a target IED of each “step range” in the inner loop is updated according to equation (13):

$\begin{matrix} {S_{l}^{*{({{h + 1},k})}} = {S_{l}^{*{({h,k})}}\frac{C^{({h,k})}}{C^{*{(k)}}}}} & (13) \end{matrix}$

where, C*^((k)) denotes a target cost updated in the kth outer loop, C^((h,k)) denotes a current cost updated in the kth outer loop and the hth inner loop. An initial target IED S_(i)*^((0,k)) of each “step range” when the process enters the inner loop is calculated by equation (14):

$\begin{matrix} {S_{i}^{*{({0,k})}} = \left\{ \begin{matrix} {{\overset{\_}{S}}_{i}^{(k)},{{\overset{\_}{S}}_{i}^{(k)} \geq {{\overset{\_}{S}}_{}^{(k)}*V_{threshold}}}} \\ {{\overset{\_}{S}}_{}^{(k)},{{\overset{\_}{S}}_{i}^{(k)} < {{\overset{\_}{S}}_{}^{(k)}*V_{threshold}}}} \end{matrix} \right.} & (14) \end{matrix}$

Where, V_(threshold) is a target IED threshold coefficient in the “step range”, S ^((k)) is an average IED of all cells in the kth outer loop; S _(i) ^((k)) is an average IED of all cells in the “step range” [id_(i−1),id_(i)] (i=1, . . . , m=1), obtained by equation (15):

$\begin{matrix} {{\overset{\_}{S}}_{i}^{(k)} = {\frac{1}{{id}_{i + 1} - {id}_{i} + 1}{\sum\limits_{{i'} = {id}_{l}}^{{id}_{i + 1}}S_{i'}^{(k)}}}} & (15) \end{matrix}$

A schematic diagram of the step target IED function constructed by the above steps is shown in FIGS. 4A-4B. A step target IED function of H_(threshold)=0 and V_(threshold)=1 is shown. in FIG. 4A, and a step target IED function of H_(threshold)=3 and V_(threshold)=1 is shown in FIG. 4B.

-   -   S4.2: Cell material updating rule based on a PID control         strategy

A cell material updating rule with a certain control strategy is to make a current cost in the inner loop converged to a target cost. The larger the flow stress of the cell material, the more difficult thin-walled frame structure is to deform in the local region, and the smaller its IED. Conversely, the smaller the flow stress of the cell material, the more easily thin-walled frame structure is to deform in the local region, and the larger its IED. A current IED of each cell is compared with the value of the SIED* function to make the current cost of the inner loop converge to the target cost: if the cell IED is lower than the SIED*, the cell material should be changed to the material with a lower flow stress.

The energy absorption capacity of thin-walled frame structures is dependent of the geometrical characteristics and material properties, in which, the key indicators affecting material properties include yield strength, tensile strength, hardening index and so on. The flow stress calculated by equation (16) can be generally employed to measure the overall material strength, which is adopted as a basis to select material.

$\begin{matrix} {\sigma_{f} = \sqrt{\frac{\sigma_{y}\sigma_{u}}{1 + n}}} & (16) \end{matrix}$

where, σ_(y) is a yield strength, σ_(u) is a tensile strength, and n=0.1 is a hardening index.

Since body material is a discrete variable, its optimization design belongs to the optimization problem with discrete variable. In addition, a specified material normally has an assured combination of different material parameters. Therefore, large amount of complex relationships among material parameters would also be introduced, which would no doubt lead to a high computational complexity of the optimization problems. To handle the difficulties mentioned above, A so-called nominal flow stress (continuous variable) is defined and updated by equations (17)-(20), which is compared with the actual flow stress of the candidate material in turn. Then the candidate material, of which the actual flow stress is closest to the nominal flow stress, is selected as the material of current cell. Finally, the material parameters of the current cell are replaced by the mechanical parameters of the selected material, i.e., density, elastic modulus, yield stress and so on.

The specific steps of cell material update are listed as follows:

-   -   S4.2.1: Definition of candidate material library and nominal         flow stress

A candidate material library of l (l≥2) materials is defined as follows:

$\begin{matrix} {{{Mat} = \left\{ {{{Mat}(1)},L,{{Mat}(s)},L,{{Mat}(l)}} \right\}},{1 \leq s \leq l}} \\ {= \left\{ {\left( {\rho_{1},E_{1},\sigma_{yi},\sigma_{u1},\sigma_{f1},L} \right),L,{\left( {\rho_{s},E_{s},\sigma_{ys},\sigma_{us},\sigma_{fs},L} \right)L},\text{ }\left( {\rho_{l},E_{l},\sigma_{yl},\sigma_{ul},{\sigma_{fl}L}} \right)} \right\}} \end{matrix}$

Where,

$\sigma_{fx} = \sqrt{\frac{\sigma_{ys}\sigma_{us}}{1 + n}}$

is a flow stress of the sth material in the candidate library, ρ_(s) is a material density of the sth material in the candidate library, and E_(s) is an elastic modulus of the sth material in the candidate library.

To solve the discrete optimal problems of body materials, we define a nominal flow stress, which is a positive non-physical parameter.

-   -   S4.2.2: The nominal flow stress for current cells is updated by         equation (17):

$\begin{matrix} {\sigma_{\Omega_{i,j}}^{\prime({{h + 1},k})} = \left\{ \begin{matrix} {\sigma_{\Omega_{i,j}}^{\min},} & {{\sigma_{\Omega_{i,j}}^{\prime({h,k})} + {\Delta\sigma_{\Omega_{i,j}}^{\prime({h,k})}}} < \sigma_{\Omega_{i,j}}^{\min}} \\ {{\text{?} + {\Delta\sigma_{\Omega_{l,j}}^{\prime({h,k})}}},} & {\sigma_{\Omega_{l,j}}^{\min} \leq {\sigma_{\Omega_{l,j}}^{\prime({h,k})} + {\Delta\sigma_{\Omega_{i,j}}^{\prime({h,k})}}} \leq \sigma_{\Omega_{i,j}}^{\max}} \\ {\sigma_{\Omega_{i,j}}^{\max},} & {{\sigma_{\Omega_{i,j}}^{\prime({h,k})} + {\Delta\sigma_{\Omega_{l,j}}^{\prime({h,k})}}} > \sigma_{\Omega_{i,j}}^{\max}} \end{matrix} \right.} & (17) \end{matrix}$ ?indicates text missing or illegible when filed

in which, σ_(Ω) _(i,j) ′^((h,k)) is a nominal flow stress of the jth cell in the hth inner loop and the kth outer loop, σ_(Ω) _(i,j) ′^((h+1,k)) is a nominal flow stress of the jth cell in the (h+1)th inner loop and the kth outer loop, σ_(Ω) _(i,j) ^(min) and σ_(Ω) _(i,j) ^(max) are respectively a minimize and a maximum value of the actual flow stress of the jth cell in the subdomain Ω_(i), Δσ_(Ω) _(i,j) ′^((h,k)) is a nominal flow stress variation of the jth cell in the subdomain Ω_(i) in the hth inner loop and the kth outer loop:

Δσ_(Ω) _(i,j) ′^((h,k))=(σ_(Ω) _(i,j) ^(max)−σ_(Ω) _(i,j) ^(min))f(e _(Ω) _(i,j) ^((h,k)))  (18)

where e_(Ω) _(i,j) ^((h,k)) denotes a difference between the current IED S_(Ω) _(i,j) ^((k)) and the target IED S_(m′)*^((h,k)) as follows

$\begin{matrix} {e_{\Omega_{i,j}}^{({h,k})} = \frac{S_{\Omega_{l,j}}^{(k)} - S_{m^{\prime}}^{*{({h,k})}}}{S_{m^{\prime}}^{*{({h,k})}}}} & (19) \end{matrix}$

f(e_(Ω) _(i,j) ^((h,k))) is a PID control function for updating the nominal flow stress, which integrates the advantages of the proportional control, integral control and derivative control to improve the control performance:

$\begin{matrix} {{f\left( e_{\Omega_{l,j}}^{({h,k})} \right)} = {{K_{p}e_{\Omega_{i,j}}^{({h,k})}} + {K_{i}\left\lbrack {\text{?} + {\sum\limits_{\tau = 1}^{k - 1}\text{?}}} \right\rbrack} + {K_{d}\left\lbrack {e_{\Omega_{i,j}}^{({h,k})} - e_{\Omega_{l,j}}^{({k - 1})}} \right\rbrack}}} & (20) \end{matrix}$ ?indicates text missing or illegible when filed

where, K_(p) is a proportional control coefficient, K_(i) is an integral control coefficient, K_(d) is a differential control coefficient, e_(Ω) _(i,j) ^((r)) is a relative deviation item of the rth outer loop, e_(Ωdi i,j) ^((k−1)) is a relative deviation item of the (k−1)th outer loop.

-   -   S4.2.3: The nominal flow stress is compared with the actual flow         stress of each material in the candidate material library and         the candidate material of which the actual flow stress is         closest to the nominal flow stress is selected as the material         of the current cell (“selected material” for short):

$\begin{matrix}  & (21) \end{matrix}$ $\left\{ \begin{matrix} {{\sigma_{\Omega_{i,j}}^{\prime({{h + 1},k})} - \sigma_{fp}} = {\min\left( {{\sigma_{\Omega_{l,j}}^{\prime({{h + 1},k})} - \sigma_{f1}},L,{\sigma_{\Omega_{i,j}}^{\prime({{h + 1},k})} - \sigma_{fs}},L,{\sigma_{\Omega_{i,j}}^{\prime({{h + 1},k})} - \sigma_{fl}}} \right)}} \\ {\sigma_{\Omega_{i,j}}^{({{h + 1},k})} = \sigma_{fp}} \\ {{Mat}_{\Omega_{i,j}}^{({{h + 1},k})} = {{Mat}(p)}} \end{matrix} \right.$

in which, p denotes a position of the selected material in the candidate material library, σ_(fp) is an actual flow stress of the selected material, σ_(Ω) _(i,j) ^((h+1,k)) is an actual flow stress of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop, Mat(p) denotes a selected material. Mat_(Ω) _(i,j) ^((h+1,k)) denotes a selected material of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop.

-   -   S4.2.4: Cell material properties are replaced by the selected         material properties

Material properties of a cell are replaced by the selected material properties mentioned above by equation (22):

$\begin{matrix} {\begin{pmatrix} \rho_{\Omega_{i,j}}^{({{h + 1},k})} \\ \text{?} \\ \sigma_{y,\Omega_{i,j}}^{({{h + 1},k})} \\ \sigma_{u,\Omega_{i,j}}^{({{h + 1},k})} \\ \text{?} \\ L \end{pmatrix} = \begin{pmatrix} \rho_{p} \\ E_{p} \\ \sigma_{yp} \\ \sigma_{up} \\ \sigma_{fp} \\ L \end{pmatrix}} & (22) \end{matrix}$ ?indicates text missing or illegible when filed

where, ρ_(Ω) _(i,j) ^((h+1,k)) is a material density jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop, E_(Ω) _(i,j) ^((h+1,k)) is an elastic modulus of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop, σ_(y,Ωi,j) ^((h+1,k)) is a yield stress of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop, σ_(u,Ωi,j) ^((h+1,k)) is a tensile strength of the jth cell in the subdomain Ω_(i) in the (h+1) inner loop and the kth outer loop, and (σ_(f,Ω) _(i,j) ^((h+1,k)) is a flow stress of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop; ρ_(p) is a material density of the selected material, E_(p) is an elastic modulus of the selected material, σ_(yp) is a yield stress of the selected material, σ_(up) is a tensile strength of the selected material, and σ_(fp) is a flow stress of the selected material.

According to the above steps, a cell cost C_(Ω) _(i,j) ^((h+1,k)) can be calculated as follows:

C _(Ω) _(i,j) ^((h+1,k))=ξ_(Ω) _(i,j) ^((h+1,k))ρ_(Ω) _(i,j) ^((h+1,k)) t _(Ω) _(i,j) A _(Ω) _(i,j)   (23)

where, C_(Ω) _(i,j) ^((h+1,k)) is a cost of the selected material of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop, ξ_(Ω) _(i,j) ^((h+1,k)) is a price of the selected material of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop, and ρ_(Ω) _(i,j) ^((h+1,k)) is a density of the selected material of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop; t_(Ω) _(i,j) is a thickness of the jth cell in the subdomain Ω_(i) and A_(Ω) _(i,j) is an area of the jth cell in the subdomain Ω_(i).

-   -   S4.3: The process returns to S4.1 if the convergence condition         of the inner loop is not satisfied or exists the inner loop and         continue to S5.

The convergence condition of the inner loop is:

|C ^((h,k)) −C* ^((k))|<ε_(i) or k _(i) ≥k _(1max)  (24)

where,

$C^{({h,k})} = {\sum\limits_{i = 1}^{l}{\sum\limits_{j = 1}^{N}C_{\Omega_{i,j}}^{({h,k})}}}$

denotes a total cost in the kth outer loop and the hth inner loop, C*^((k)) is a target cost defined in the kth outer loop, ε₁ is a cost convergence factor, k₁ denotes the number of iterations in the inner loop, and k_(1max) denotes a maximum number of iterations in the inner loop.

-   -   S5: The optimal results are written out if the global         convergence conditions in the outer loop are satisfied or the         M-SHCA algorithm returns to S3.

The M-SFICA algorithm will be terminated if one of the following three convergence conditions is satisfied:

-   -   (1) The number of iterations k (namely, the number of finite         element simulation analyses) in the outer loop exceeds a user         defined a maximum number of iterations k_(max).     -   (2) p_(f)>p_(f)* and a current design point is a feasible         solution, in which P_(f) denotes the number of iterations where         infeasible solutions continuously appear and p_(f)* denotes a         maximum number of iterations where infeasible solutions         continuously appear.     -   (3) The difference of design variables between two iterations is         very small, namely:

${\sum\limits_{i = 1}^{N}{❘{\text{?} - \text{?}}❘}} < \varepsilon_{2}$ ?indicates text missing or illegible when filed

in which, N is the total number of cells, and ε₂ represents a global convergence factor.

Embodiment

To verify the convergence and efficiency of the M-SHCA algorithm, it is employed to optimize the material distribution and cost of a car body frame under side collisions. The total weight of the full-vehicle crash FE model is 1346 kg including 276838 elements and 284961 nodes, in which body in white (BIW) adopts the shell elements and engine, gearbox, suspension system, etc. adopt the solid elements. The piecewise elastoplastic materials are used for the deformable structures and the rigid materials are adopted for the undeformable structures. The automatic single surface, automatic surface to surface, automatic node to surface algorithms are defined for the possible contact positions during side collisions. According to the requirements of the regulation titled “The protection of the occupants in the event of a lateral collision” (GB 20071-2006), a mobile deformable barrier (MDB) with a weight of 950 kg should hit a target vehicle perpendicularly at an initial velocity of 50 km/h, as shown in FIG. 5 .

Step 1: Definition of the Subdomain CA Model and Design Variable

During the vehicle side collisions, B-pillar, sill, doors, and roof middle crossbeam appear large deformation which are the main energy absorbing structures and A-pillar, roof rail, seat crossbeam, and roof crossbeam are the main structures to transfer impact loading. Therefore, the material of 34 parts of 14 assemblies, such as the A-pillar, B-pillar, sill, roof rail, front and rear doors, rear side member, seat crossbeam, and roof crossbeam are defined as the design variables,

The detailed steps to define the subdomain. CA model for the car body beam frame are provided as follows:

-   -   Step (1): Subdomain fragmentation: the design space is         fragmented into several independent subdornains denoted as Ω₁         based on the topological connection characteristics of the car         body frame structure. For example, the assemblies such as         A-pillar, B-pillar, and sill, are respectively defined as a         subdomain Ω₁, a subdomain Ω₂, and a subdomain Ω₃, as shown in         FIG. 6 .     -   Step (2): Cell definition: Each component is defined as a cell         Ω_(i,j) in the subdomain Ω₁(i=1,2,L,l), where the subscript i of         Ω_(i,j) denotes an index of the ith subdomain, and the subscript         j denotes a location of current cell in the ith subdomain. The         subscripts j of each cell Ω_(i,j) are sequentially numbered from         small to large according to principles of from inside to         outside, from front to back, and from bottom to top.     -   Step (3): Cell state variable definition: The design variables         (e.g. material) and field variables (e.g. IEDs) are sequentially         defined for each cell.     -   Step (4). cell neighborhood definition: All subdomains are         traversed, for cells in the same subdomain, the neighboring         cells of the current cell are determined according to orders of         subscripts j, and a set of the neighboring cells of the current         cell is referred to as the cell neighborhood. Three subdomains         are defined in FIG. 6 . All cells within the cell radius of r=1         from the center of the current cell are referred to as the         neighboring cells of the current cell. The number of neighboring         cells of a current cell Ω_(1,2) is 2, ({circumflex over (N)}_(Ω)         _(1,2) =2) the number of neighboring cells of a current cell         Ω_(2,1) is 1 ({circumflex over (N)}_(Ω) _(2,1) =1), and the         number of neighboring cells of a current cell Ω_(3,1) is 1         ({circumflex over (N)}_(Ω) _(3,1) =1).

Following the above 4 steps, a total of 14 subdomains and a total of 34 thickness variables are defined for the car body beam frame model, as shown in TABLE 2. The material parameters of the candidate material library are listed in FIG. 7 and TABLE 3.

TABLE 2 The SCA model and the design variables of body beam frame SCA CA Design variables Name Sy Name Symbol Symbol Initial A-pillar Ω₁ A-pillar inner panel Ω_(1, 1) Mat₁ HC340 A-pillar reinforcement#1 Ω_(1, 2) Mat₂ B280VK A-pillar reinforcement#2 Ω_(1, 3) Mat₃ HC340 A-pillar outer panel Ω_(1, 4) Mat₄ HC420 B-pillar Ω₂ B-pillar inner panel Ω_(2, 1) Mat₅ DP980 B-pillar reinforcement Ω_(2, 2) Mat₆ DP980 Sill Ω₃ Sill inner panel#1 Ω_(3, 1) Mat₇ HC420 Sill inner panel#2 Ω_(3, 2) Mat₈ HC420 Sill reinforcement Ω_(3, 3) Mat₉ HC420 Roof rail Ω₄ A-pillar roof rail#1 Ω_(4, 1) Mat₁₀ B210P1 A-pillar roof rail#2 Ω_(4, 2) Mat₁₁ HC420 B-pillar roof rail Ω_(4, 3) Mat₁₂ HC420 C-pillar roof rail Ω_(4, 4) Mat₁₃ HC420 Front door Ω₅ Front door anti-collision Ω_(5, 1) Mat₁₄ HC340 Front door anti-collision Ω_(5, 2) Mat₁₅ HC420 Front door anti-collision Ω_(5, 3) Mat₁₆ HC340 Front door inner panel Ω_(5, 4) Mat₁₇ HC420 Rear door Ω₆ Rear door anti-collision beam Ω_(6, 1) Mat₁₈ HC340 Rear door anti-collision Ω_(6, 2) Mat₁₉ HC420 Rear door anti-collision beam Ω_(6, 3) Mat₂₀ HC340 Rear door inner panel Ω_(6, 4) Mat₂₁ HC420 Rear door anti-collision beam Ω_(6, 5) Mat₂₂ HC340 Rear door anti-collision Ω_(6, 6) Mat₂₃ HC420 Rear door anti-collision beam Ω_(6, 7) Mat₂₄ HC340 mounting panel#4 Rear side member Ω₇ Rear side member inner panel Ω_(7, 1) Mat₂₅ HC420 Rear side member outer Ω_(7, 2) Mat₂₆ HC420 Seat crossbeam Ω₈ Seat crossbeam lining panel Ω_(8, 1) Mat₂₇ HC420 Seat crossbeam Ω_(8, 2) Mat₂₈ HC420 Front side member Ω₉ Front side member rear Ω_(9, 1) Mat₂₉ HC420 Seat rear crossbeam Ω₁₀ Seat rear crossbeam Ω_(10, 1) Mat₃₀ HC420 Rear floor front Ω₁₁ Rear floor front crossbeam Ω_(11, 1) Mat₃₁ HC420 Roof front Ω₁₂ Roof front crossbeam Ω_(12, 1) Mat₃₂ HC420 Roof middle Ω₁₃ Roof middle crossbeam Ω_(13, 1) Mat₃₃ DP980 Roof rear crossbeam Ω₁₄ Roof rear crossbeam Ω_(14, 1) Mat₃₄ HC420

TABLE 3 Material parameters for the candidate material library Density Elastic Yield Tensile Hardening Flow Price Material (g/cm³) modulus strength strength Index stress (CNY/kg) DC01 7.8 210 209 386 0.1 270.8 7 B170P1 7.8 210 221 470 0.1 307.3 8 B210P1 7.8 210 280 519 0.1 363.5 9 B280VK 7.8 210 363 632 0.1 456.7 10 HC340 7.8 210 398 809 0.1 541.0 11 HC420 7.8 210 550 935 0.1 683.7 12 980DP 7.8 210 720 1167 0.1 874.0 13 1180DP 7.8 210 1009 1167 0.1 1112.8 14 B1500HS 7.8 210 1122 1667 0.1 1304.0 15

Step 2: Definition of Output Response

In the side collisions simulation of the car body frame, B-pillar is a key component that resists excessive deformation of the body structure and reduces the speed of body intrusion, Excessive deformation of B-pillar will lead the body structure to invade a passenger compartment significantly while reducing the living space of the passenger compartment and causing crash injuries to the passenger. The soft -tissue organs such as heart and lungs of passenger are very sensitive to speed changing of the chest position, If the intrusion velocity is too high, the vital organs in chest will be damaged seriously. Therefore, the maximum intrusion amounts and maximum intrusion velocities of B-pillar corresponding to the chest and pelvic positions are respectively selected as the crashworthiness indexes and output responses under side collisions, which are denoted as d₁(Mat), v₁(Mat), d₂(Mat) and v₂(Mat), respectively. As depicted in FIG. 8 , a measuring point of B-pillar corresponding to chest is the point A on the inner panel of B-pillar parallel to passenger's chest and a measuring point of B-pillar corresponding to pelvis is the point B on the inner panel of B-pillar parallel to passenger's pelvis.

Step 3: Definition of Optimization Formulation

In this embodiment, the total cost of 34 parts in TABLE 2 is used as the objective function in which the initial cost is CNY 1291. The maximum invasion amounts and maximum intrusion speeds corresponding to the measuring points at B-pillar (points A and B) are defined as the constraint functions, in which the maximum intrusions at points A and B are 204.70 mm of 270.30 mm, respectively; the maximum initial intrusion velocities at points A and B are 7.78 m/s and 7.85 m/s, respectively. To make the initial full vehicle model meet the requirements of GB 20071-2006, the maximum intrusion amount and the maximum intrusion velocity should be less than or equal to 180 mm and 7.50 m/s, respectively. The initial values and design goals of the output response of B-pillar corresponding to the chest and pelvic positions are shown in TABLE 4, and the corresponding optimization equation is given as follows:

$\begin{matrix} {\left\{ \begin{matrix} \min & {{{Cos}{t({Mat})}} = {\text{?}\xi_{i}\rho_{i}t_{i}A_{i}}} \\ {s.t.} & {{d_{1}({Mat})} \leq 180} \\  & {{d_{2}({Mat})} \leq 180} \\  & {{v_{1}({Mat})} \leq 7.5} \\  & {{v_{2}({Mat})} \leq 7.5} \\  & {{Mat}_{i} \in \begin{pmatrix} {{{DC}01},{B170P1},{B210P1},{B280{VK}},{{HC}340},} \\ {{{HC}420},{980{DP}},{1180{DP}},{B1500{HS}}} \end{pmatrix}} \\  & {{Mat} = \left( {\text{?},\ldots,{Mat}_{i},\ldots,{Mat}_{34}} \right)} \end{matrix} \right.} & (25) \end{matrix}$ ?indicates text missing or illegible when filed

in which, ξ_(i) is price of the ith cell material, ρ_(i) is material density of the ith cell material, t_(i) is thickness of the ith cell, A_(i) is area of the ith cell, Mat_(i) is material number of the ith cell, and DC01, B170P1, B210P1, B280VK, HC340, HC420, 980DP, 1180DP, B1500HS are candidate materials in the material library.

TABLE 4 Output responses and design target Responses Symbol Baseline design Design target Total cost of the 34 parts (CNY) Cost(Mat) 1291 minimize B-pillar maximum intrusion amount at chest d₁(Mat) 207.30 ≤180 location (mm) B-pillar maximum intrusion amount at pelvis d₂(Mat) 204.70 ≤180 location (mm) B-pillar maximum intrusion velocity at chest v₁(Mat) 7.78 ≤7.5 location (m/s) B-pillar maximum intrusion velocity at pelvis v₂(Mat) 7.85 ≤7.5 location (m/s)

Step 4: Optimization Results and Discussion

The horizontal IED target (HIED*) function is commonly calculated by the conventional HCA method in each iteration of the inner loop to make the current cost of the inner loop converged to the target mass by updating material with the PID control strategy.

In this embodiment, the M-SHCA algorithm adopting the. HIED* function for cell material updating in the inner loop is referred to as “M-SHCA#1” and the M-SHCA algorithm adopting the SIED* function for cell material updating in the inner loop is referred to as “M-SHCA#2”. To validate the convergence and efficiency of T-SHCA#2, the optimization equation in Equation (25) is separately solved by T-SHCA#1, T-SHCA#2 and parallel EGO-PCEI. The optimization results and the FEAs' numbers of the three algorithms are compared. The detail parameters used by T-SHCA#1, T-SHCA#2 are listed in TABLE 5, while those of parallel EGO-PCEI are listed in TABLE 6.

TABLE 5 Parameters configuration for M-SHCA#1 and M-SHCA#2 Parameter name Symbol M-SHCA#1 M-SHCA#2 Total number of cells N 34 34 Number of subdomains l 14 14 Cell radius r 1 1 Width threshold of “step interval” H_(threshold) — 4 IED target threshold coefficient of “step interval” V_(threshold) — 1.1 Proportional control coefficient of cell material K_(p) 0.03 0.03 variation Integral control coefficient of cell material variation K_(i) 0.0001 0.0001 Differential control coefficient of cell material K_(d) 0.0001 0.0001 variation Scale factor of the penalty value of the target cost K_(q) 0.15 0.15 Maximum penalty of target cost (CNY) ΔC 35 35 Maximum iteration number for the consecutive p*_(f) 10 10 infeasible solutions Maximum number of iterations in inner loop k_(1max) 2000 2000 Maximum number of iterations in outer loop k_(max) 50 50 Inner loop convergence factor ε₁ 0.001 0.001 Outer loop convergence factor ε₂ 0.001 0.001

TABLE 6 Parameters configuration for parallel EGO-PCEI Parameter name Symbol Parameter value Number of design variables N_(v) 34 Number of initial samples N_(initial) 4 Number of parallel calculations N_(parallel) 4 Maximum number of FEAs Max_(FEA) 300

The iteration processes of M-SHCA#1 are illustrated in FIGS. 9A-9C, in which M-SHCA#1 converges after 23 iterations, an optimal solution is found in the 12th iteration, and a total of 23 FEAs are conducted. The iteration processes of M-SHCA#2 are illustrated in FIGS. 10A-10B, in which M-SHCA#2 converges after 18 iterations, an optimal solution is found in the 8th iteration, and a total of 18 FEAs are conducted. With the comparisons between FIGS. 9A-9C and FIGS. 10A-10B, M-SHCA#1 cannot find the global optimal solution, and it is very hard to decrease the total cost of body frame; M-SHCA#2 converges when its objective function quickly drops to CNY1101. The iteration processes of parallel EGO-PCEI are illustrated in FIGS. 11A-11B, in which the global optimal solution is found after 67 iterations and a total of 300 FEAs are conducted before converging. Table 7 lists the initial designs, the responses of the optimal solutions and the total number of FEAs when the responses of T-SHCA#1 , T-SHCA#2, and parallel EGO-PCEI converge to the optimal solutions. The computational resources and computational time spent by the three methods are compared in FIG. 12 .

TABLE 7 Comparison between the initial design and the optimal solutions Initial M- M- Parallel EGO- design SHCA#1 SHCA#2 PCEI Cost(Mat) (CNY) 1291 1142 1101 1051 d₁(Mat) (mm) 207.30 168.3 176.20 177.20 d₂(Mat) (mm) 204.70 179.4 180.00 177.50 v₁(Mat) (m/s) 7.78 6.97 6.96 6.87 v₂(Mat) (m/s) 7.85 7.33 7.28 6.74 Number of FEAs — 12 8 271 when converging to optimal solution Total number of 1 23 18 300 FEAs Number of CPUs 8 8 8 32 Total time 3 69 54 225 consumed (h)

With the comparison and analysis of FIGS. 9A-12 and TABLE 7, (1) M-SHCA#1 and M-SHCA#2 not only spend less computing resources, but also have higher optimization efficiency; (2) M-SHCA#1 is very easy to trap into a local optimal solution, but very hard to decrease the total cost of body frame; (3) M-SHCA#2 not only has a high optimization efficiency, but also a stronger capability for global searching and cost decreasing; (4) the global optimization precision of M-SHCA#2 is very close to that of parallel EGO-PCEI, but the optimization efficiency of M-SHCA#2 is obviously higher than that of parallel EGO-PCEI. In summary, M-SHCA#2 has higher global searching efficiency and numerical precision for solving the time-consuming nonlinear dynamic response optimization problems with multiple discrete variables and its effectiveness has been validated by engineering case.

The optimization effect of the M-SHCA algorithm based on the SIED* function (i.e. M-SHCA#2) is further discussed here. The total cost of 34 parts and the performance improvement percentage under side collisions before and after optimization as listed in Table 8, in which the optimal solution obtained by the M-SHCA algorithm has achieved a cost reduction effect of 14.72%, while d₁(Mat), d₂(Mat), v₁(Mat) and v₂(Mat) have reduced by 15.00%, 12.07%, 10.54%, and 7.26%, respectively. The proposed algorithm not only reduces the total cost of body frame to a large extent, but also significantly improves the safety of vehicle under side collisions.

TABLE 8 Performance improvement percentages before and after optimization Cost d₁ d₂ v₁ v₂ (Mat) (Mat) (Mat) (Mat) (Mat) (CNY) (mm) (mm) (m/s) (m/s) Initial design 1291 207.30 204.70 7.78 7.85 Optimal solution 1101 176.20 180.00 6.96 7.28 Relative change rate (%) −14.72% −15.00% −12.07% −10.54% −7.26%

The material distributions of the optimization solution and the initial body frame are compared in TABLE 9, in which the material of the optimal solution have been distributed more reasonably compared with the initial design; with the optimal material distribution for body frame, the side collision safety performance of body frame can be improved while the total cost of body frame can be greatly reduced. In other words, the multi-material body frame has a stronger potential to improve its crash safety and reduce its total cost than initial body frame with a small amount of materials.

TABLE 9 Material distribution of body frame before and after optimization Initial Optimal Name Symbol material material A-pillar inner panel Mat₁ HC340 DC01 A-pillar reinforcement#1 Mat₂ B280VK DC01 A-pillar reinforcement#2 Mat₃ HC340 DC01 A-pillar outer panel Mat₄ HC420 DC01 B-pillar inner panel Mat₅ DP980 B1500HS B-pillar reinforcement Mat₆ DP980 B1500HS Sill inner panel#1 Mat₇ HC420 1180DP Sill inner panel#2 Mat₈ HC420 B1500HS Sill reinforcement Mat₉ HC420 HC420 A-pillar roof rail#1 Mat₁₀ B210P1 DC01 A-pillar roof rail#2 Mat₁₁ HC420 DC01 B-pillar roof rail Mat₁₂ HC420 DC01 C-pillar roof rail Mat₁₃ HC420 DC01 Front door anti-collision beam mounting Mat₁₄ HC340 B1500HS panel#1 Front door anti-collision beam Mat₁₅ HC420 B1500HS Front door anti-collision beam mounting Mat₁₆ HC340 B1500HS panel#2 Front door inner panel reinforcement Mat₁₇ HC420 DC01 Rear door anti-collision beam mounting panel#1 Mat₁₈ HC340 DC01 Rear door anti-collision beam#1 Mat₁₉ HC420 DC01 Rear door anti-collision beam mounting panel#2 Mat₂₀ HC340 DC01 Rear door inner panel reinforcement Mat₂₁ HC420 DC01 Rear door anti-collision beam mounting panel#3 Mat₂₂ HC340 DC01 Rear door anti-collision beam#2 Mat₂₃ HC420 DC01 Rear door anti-collision beam mounting panel#4 Mat₂₄ HC340 DC01 Rear side member inner panel Mat₂₅ HC420 B280VK Rear side member outer panel Mat₂₆ HC420 B280VK Seat crossbeam lining panel Mat₂₇ HC420 B1500HS Seat crossbeam Mat₂₈ HC420 B1500HS Front side member rear section Mat₂₉ HC420 DC01 Seat rear crossbeam Mat₃₀ HC420 B210P1 Rear floor front crossbeam Mat₃₁ HC420 B170P1 Roof front crossbeam Mat₃₂ HC420 B170P1 Roof middle crossbeam Mat₃₃ DP980 B1500HS Roof rear crossbeam Mat₃₄ HC420 DC01

From FIGS. 13A-13B, it is concluded that (1) the intrusion amounts and the intrusion velocity of the optimal solution are greatly improved compared with the initial design, (2) B-pillar and roof middle crossbeam have been strengthened in the optimal solution, of which the deformation modes have been significantly improved compared with the initial design.

From the discussion mentioned above, it is concluded that the M-SHCA algorithm based on the SIED* function have a higher efficiency of global searching than that based on the HIED* function for solving the large scale nonlinear dynamic responses structural optimization problems with many discrete design variables, So the M-SHCA algorithm based on the SIED* function can be employed to effectively solve the optimization problems including the intrusion amounts and the intrusion velocity constraints, specially the nonlinear dynamic structural optimization problems with large scale discrete design variables.

The described embodiment is a preferred embodiment of the present invention, but the present invention is not limited to the aforementioned embodiment. Any obvious improvements, substitutions or modifications that can be made by those skilled in the art without departing from the essential content of the present invention shall fall within the protection scope of the present invention. 

What is claimed is:
 1. A material-based subdomain hybrid cellular automata algorithm for solving material optimization of thin-walled frame structures, comprising the following steps: S1: establishing an initial crash finite element model, constructing a subdomain cellular automata model, defining material variables and field variables of the thin-walled frame structures, and employing the initial crash finite element model for material and cost optimization; S2: executing an outer loop: calculating a cell internal energy density and a constraint value at a current design point by finite element analysis and updating a target cost by a penalty function method according to an extent of the current design point violating a constraint boundary; S3: executing an inner loop with the following steps: S3.1: constructing a step internal energy density (IED) target (SIED*) function and updating a target IED; S3.2: updating a cell material by a material updating rule based on a proportional integral derivative (PID) control strategy; specifically: defining a candidate material library and a nominal flow stress of each material, updating a nominal flow stress of a current cell, comparing the nominal flow stress with a true flow stress of each material in the candidate material library, and selecting a candidate material closest to the nominal flow stress as a selected material of the current cell, replacing material parameters of the current cell with mechanical parameters of the selected material; S3.3: executing S4 and exiting an inner loop if the inner loop is convergent, otherwise returning to S3.1; S4: writing out optimal results if global convergence conditions in the outer loop are satisfied, otherwise returning to S2 for updating the cell material in the inner loop.
 2. The material-based subdomain hybrid cellular automata algorithm for solving the material optimization of the thin-walled frame structures according to claim 1, wherein the candidate material library is defined as follows: $\begin{matrix} {{{Mat} = \left\{ {{{Mat}(1)},L,{{Mat}(s)},L,{{Mat}(l)}} \right\}},{1 \leq s \leq l}} \\ {= \left\{ {\left( {\rho_{1},E_{1},\sigma_{y1},\sigma_{u1},\sigma_{f1},L} \right),L,{\left( {\rho_{s},E_{s},\sigma_{ys},\sigma_{us},\sigma_{fs},L} \right)L},} \right.} \\ \left. {}\left( {\rho_{l},E_{l},\sigma_{yl},\sigma_{ul},\sigma_{fl},L} \right) \right\} \end{matrix}$ wherein, ρ_(s) is a density of a sth material in the candidate material library; E_(s) is an elastic modulus of the sth material in the candidate material library; σ_(ys) is a yield strength of the sth material in the candidate material library; $\sigma_{fs} = \sqrt{\frac{\sigma_{ys}\sigma_{us}}{1 + n}}$ is a flow stress of the sth material in the candidate material library; σ_(us) is a tensile strength of the sth material in the candidate material library; l≥2 is a number of materials in the candidate library.
 3. The material-based subdomain hybrid cellular automata algorithm for solving the material optimization of the thin-walled frame structures according to claim 1, wherein the nominal flow stress is a non-physical parameter, and the non-physical parameter is a positive real number.
 4. The material-based subdomain hybrid cellular automata algorithm for solving the material optimization of the thin-walled frame structures according to claim 1, wherein the nominal flow stress of each cell is updated by a following equation: $\sigma_{\Omega_{i,j}}^{\prime({{h + 1},k})} = \left\{ \begin{matrix} {\sigma_{\Omega_{i,j}}^{\min},} & {{\sigma_{\Omega_{i,j}}^{\prime({h,k})} + {\Delta\sigma_{\Omega_{i,j}}^{\prime({h,k})}}} < \sigma_{\Omega_{i,j}}^{\min}} \\ {{\sigma_{\Omega_{i,j}}^{\prime({h,k})} + {\Delta\sigma_{\Omega_{i,j}}^{\prime({h,k})}}},} & {\sigma_{\Omega_{i,j}}^{\min} \leq {\sigma_{\Omega_{i,j}}^{\prime({h,k})} + {\Delta\sigma_{\Omega_{i,j}}^{\prime({h,k})}}} \leq \sigma_{\Omega_{i,j}}^{\max}} \\ {\sigma_{\Omega_{i,j}}^{\max},} & {{\text{?} + {\Delta\sigma_{\Omega_{i,j}}^{\prime({h,k})}}} > \text{?}} \end{matrix} \right.$ ?indicates text missing or illegible when filed wherein, σ_(Ω) _(i,j) ′^((h,k)) is a nominal flow stress of a jth cell in a subdomain Ω₁ in a hth inner loop and a kth outer loop; σ_(Ω) _(i,j) ′^((h+1,k)) is a nominal flow stress of the jth cell in the subdomain Ω₁ in a (h+1)th inner loop and the kth outer loop; σ_(Ω) _(i,j) ^(min) and σ_(Ω) _(i,j) ^(max) are a minimize value and a maximum value of an actual flow stress of the jth cell in the subdomain Ω₁, respectively; Δσ_(Ω) _(i,j) ′^((h,k)) is a nominal flow stress variation of the jth cell in the subdomain Ω₁ in the hth inner loop and the kth outer loop: Δσ_(Ω) _(i,j) ′^((h,k))=(σ_(Ω) _(i,j) ^(max)−σ_(Ω) _(i,j) ^(min))f(e _(Ω) _(i,j) ^((h,k))) wherein, e_(Ω) _(i,j) ^((h,k)) denotes a difference between a current IED S_(Ω) _(i,j) ^((k)) and a target IED S_(m′)*^((h,k)) and a PID control function for updating the nominal flow stress f(e_(Ω) _(i,j) ^((h,k))) is given as follows: ${f\left( e_{\Omega_{i,j}}^{({h,k})} \right)} = {{K_{p}e_{\Omega_{i,j}}^{({h,k})}} + {K_{i}\left\lbrack {e_{\Omega_{i,j}}^{({h,k})} + {\sum\limits_{\tau = 1}^{k - 1}e_{\Omega_{i,j}}^{(\tau)}}} \right\rbrack} + {K_{d}\left\lbrack {e_{\Omega_{i,j}}^{({h,k})} - e_{\Omega_{i,j}}^{({k - 1})}} \right\rbrack}}$ wherein, K_(p) is a proportional control coefficient; K_(i) is an integral control coefficient; K_(d) is a differential control coefficient; e_(Ω) _(i,j) ^((r)) is a relative deviation item of a rth outer loop; e_(Ω) _(i,j) ^((k−1)) is a relative deviation item of a (k-1)th outer loop.
 5. The material-based subdomain hybrid cellular automata algorithm for solving the material optimization of the thin-walled frame structures according to claim 4, wherein the following equations are used to select the candidate material closest to the nominal flow stress as the selected material of the current cell by comparing the nominal flow stress with the true flow stress of each material in the candidate material library: $\left\{ \begin{matrix} {{\sigma_{\Omega_{i,j}}^{\prime({{h + 1},k})} - \sigma_{fp}} = {\min\left( {{\sigma_{\Omega_{i,j}}^{\prime({{h + 1},k})} - \sigma_{f1}},L,{\sigma_{\Omega_{i,j}}^{\prime({{h + 1},k})} - \sigma_{fs}},L,{\sigma_{\Omega_{i,j}}^{\prime({{h + 1},k})} - \sigma_{fl}}} \right)}} \\ {\text{?} = \sigma_{fp}} \\ {{Mat}_{\Omega_{i,j}}^{({{h + 1},k})} = {{Mat}(p)}} \end{matrix} \right.$ ?indicates text missing or illegible when filed wherein, p denotes a position of the selected material in the candidate material library; σ_(fp) is an actual flow stress of the selected material; σ₁₀₆ _(i,j) ^((h+1,k)) is an actual flow stress of the jth cell in the subdomain Ω₁ in the (h+1)th inner loop and the kth outer loop; Mat(p) denotes the selected material; Mat_(Ω) _(i,j) ^((h+1,k)) denotes a selected material of the jth cell in the subdomain Ω₁ in the (h+1)th inner loop and the kth outer loop; σ_(s) is a flow stress of a sth material in the candidate material library.
 6. The material-based subdomain hybrid cellular automata algorithm for solving the material optimization of the thin-walled frame structures according to claim 1, wherein the following equation is employed to replace material parameters of the current cell with mechanical parameters of the selected material: $\begin{pmatrix} \text{?} \\ E_{\Omega_{i,j}}^{({{h + 1},k})} \\ \text{?} \\ \sigma_{u,\Omega_{i,j}}^{({{h + 1},k})} \\ \sigma_{f,\Omega_{i,j}}^{({{h + 1},k})} \\ L \end{pmatrix} = \begin{pmatrix} \rho_{p} \\ E_{p} \\ \sigma_{yp} \\ \sigma_{up} \\ \sigma_{fp} \\ L \end{pmatrix}$ ?indicates text missing or illegible when filed wherein, ρ_(Ω) _(i,j) ^((h+1,k)) is a material density of a jth cell in a subdomain Ω_(i) in a (h+1)th inner loop and a kth outer loop; E_(Ω) _(i,j) ^((h+1,k)) is an elastic modulus of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop; σ_(y,Ω) _(i,j) ^((h+1,k)) is a yield tress of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop; σ_(u,Ωdi i,j) ^((h+1,k)) is a tensile strength of the jth cell in the subdomain Ω_(i) in the (h+1)th inner loop and the kth outer loop; σ_(f,Ω) _(i,j) ^((h+1,k)) is a flow stress of the jth cell in the subdomain Ω_(i) in the (h+1 )th inner loop and the kth outer loop; ρ_(p) is a material density of the selected material; E_(p) is an elastic modulus of the selected material; σ_(yp) is a yield tress of the selected material; σ_(up) is a tensile strength of the selected material; σ_(fp) is a flow stress of the selected material.
 7. The material-based subdomain hybrid cellular automata algorithm for solving the material optimization of the thin-walled frame structures according to claim 1, wherein the global convergence conditions comprise: |C ^((h,k)) −C* ^((k))|<ε₁ or k ₁ ≥k _(1max) wherein, $C^{({h,k})} = {\sum\limits_{i = 1}^{l}{\sum\limits_{j = 1}^{N}C_{\Omega_{i,j}}^{({h,k})}}}$  is a total cost in a hth inner loop and a kth outer loop, here C_(Ω) _(i,j) ^((h,k)) is a cell cost in a inner loop and a outer loop; C*^((k)) is a target cost defined in the kth outer loop; ε₁ is an inner loop convergence factor; k₁ is a number of iterations in the inner loop; k_(1max) is a maximum number of iterations in the inner loop; a cell cost of the selected material of a jth cell in a subdomain Ω₁ in a (h+1)th inner loop and the kth outer loop C_(Ω) _(i,j) ^((h+1,k)) is calculated as follows: C _(Ω) _(i,j) ^((h+1,k))=ξ_(Ω) _(i,j) ^((h+1,k))ρ_(Ω) _(i,j) ^((h+1,k)) t ₁₀₆ _(i,j) A _(Ω) _(i,j) wherein, ξ_(Ω) _(i,j) ^((h+1,k)) is a price of the selected material of the.jth cell in the subdomain Ω₁ in the (h+1)th inner loop and the kth outer loop, ρ_(Ω) _(i,j) ^((h+1,k)) is a density of the selected material of the jth cell in the subdomain Ω₁ in the (h+1)th inner loop and the kth outer loop, t_(Ω) _(i,j) is a thickness of the jth cell in the subdomain Ω₁, and A_(Ω) _(i,j) is an area of the jth cell in the subdomain Ω₁. 